##########################################
##########################################
##########################################
#####
#####              Chase Joyner
#####             882 Homework 2
#####           September 28, 2016
#####
##########################################
##########################################
##########################################

##########################
##### Problem 2 Code #####

## Part b ##
mu1 = 5
mu2 = -2.7
sigma1sq = 0.8
sigma2sq = 0.23
n1 = n2 = 10000
X1 = rnorm(n1, mu1, sqrt(sigma1sq))
X2 = rnorm(n2, mu2, sqrt(sigma2sq))
Y = X1 + X2
hist(Y, prob = TRUE, xlab = "Y", ylab = "p(Y)")
lines(density(Y), lwd = 2)
sigmasq = sigma1sq + sigma2sq
mu = mu1 + mu2
tmp = seq(-2, 6, by = 0.001)
ptmp = 1 / sqrt(2*pi*sigmasq) * exp(-1/(2*sigmasq)*(tmp - mu)^2)
lines(tmp, ptmp, col = 2, lty = 2, lwd = 2)
mean(Y)
var(Y)

## Part c ##
mu1 = 5
mu2 = -2.7
sigma1sq = 0.8
sigma2sq = 0.23
y = seq(-2, 6, by = 0.001)
N = length(y)
dY = rep(-99, N)
for(i in 1:N){
	X2.samp = rnorm(N, mu2, sqrt(sigma2sq))
	dY[i] = mean(dnorm(y[i]-X2.samp, mu1, sqrt(sigma1sq)))
}
lines(y, dY, col = "blue", lty = 5)

## Part d ##
alpha1 = 2.3
beta1 = 3
alpha2 = 2.8
beta2 = 2.2
n1 = n2 = 10000
X1 = rgamma(n1, alpha1, beta1)
X2 = rgamma(n2, alpha2, beta2)
Y = X1 + X2
hist(Y, prob = TRUE, xlab = "Y", ylab = "p(Y)", ylim = c(0,0.6))
lines(density(Y), lwd = 2)
y = seq(0, 5, by = 0.001)
N = length(y)
dY = rep(-99, N)
for(i in 1:N){
	X2.samp = rgamma(N, alpha2, beta2)
	dY[i] = mean(dgamma(y[i]-X2.samp, alpha1, beta1))
}
lines(y, dY, col = "blue", lty = 2)


x1 = rgamma(10000, alpha1, beta1)
x2 = rgamma(10000, alpha2, beta2)

hist(x1, prob = TRUE)
windows()
hist(x2, probe = TRUE)

##########################
##### Problem 3 Code #####

## Part c ##
playcraps = function(){
	num = append(1:6, 5:1)
	first.roll = sample(2:12, 1, prob = num / 36)
	if(first.roll %in% c(2,3,12)){
		return(0)
	} else if(first.roll %in% c(7,11)){
		return(1)
	} else{
		point = first.roll
		next.roll = sample(2:12, 1, prob = num / 36)
		while(next.roll != point && next.roll != 7){
			next.roll = sample(2:12, 1, prob = num / 36)
		}
		if(next.roll == point){
			return(1)
		}
		return(0)
	}
}

n = 10000
res = rep(-99, n)
for(i in 1:n){
	res[i] = playcraps()
}
mean(res)	# Should be around 49%.  Checks out.


sims = 1000000	# Number of times simulated
n = 10		# Number of games to play
i.bet = 10	# Initial bet of 10 dollars
earnings = rep(0, sims)

for(t in 1:sims){
	bet = i.bet
	tmp = 0
	for(i in 1:n){
		res = playcraps()
		if(res == 0){
			tmp = tmp - bet
			bet = 2*bet
		}
		if(res == 1){
			tmp = tmp + bet
			bet = i.bet
		}
	}
	print(t)
	earnings[t] = tmp
}
mean(earnings)
hist(earnings)








##########################
##### Problem 4 Code #####

## Part c ##
elk = function(x,p){
	val = choose(29,x[1])*p[1]^x[1]*choose(29-x[1],x[2])*p[2]^x[2]*choose(29-x[1]-x[2],x[3])*p[3]^x[3]*p[4]^x[4]
	return(val)
}

X = c(5,11,6,7)
LpA = elk(X,c(1/3,1/3,1/6,1/6))
LpB = elk(X,c(1/6,1/3,1/3,1/6))
LpC = elk(X,c(1/6,1/6,1/6,1/2))

test.stat = LpC / max(LpA, LpB, LpC)

sims = 10000
stats = rep(-99, sims)
pA = c(1/3,1/3,1/6,1/6)
pB = c(1/6,1/3,1/3,1/6)
pC = c(1/6,1/6,1/6,1/2)
for(i in 1:sims){
	rolls = sample(1:4, 29, replace = TRUE, prob = pC)
	ones = length(which(rolls == 1))
	twos = length(which(rolls == 2))
	threes = length(which(rolls == 3))
	fours = length(which(rolls == 4))
	X = c(ones,twos,threes,fours)
	mle = max(elk(X,pA),elk(X,pB),elk(X,pC))
	stats[i] = elk(X,pC) / mle
}
mean(stats <= test.stat)